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Abstract 

Statistical data is often analyzed as a contingency table, sometimes 
with empty cells called zeros. Such sparse tables can be due to scarse 
observations classified in numerous categories, as for example in ge- 
netic association studies. Thus, classical independence tests involving 
Pearson's chi-square statistic Q or Kullback's minimum discrimination 
information statistic G cannot be applied because some of the expected 
frequencies are too small. More generally, we consider goodness of fit 
tests with composite hypotheses for sparse multinomial vectors and 
suggest simple corrections for Q and G that improve and generalize 
known procedures such as Ku's. We show that the corrected statistics 
share the same asymptotic distribution as the initial statistics. We 
produce Monte Carlo estimations for the type I and type II errors on a 
toy example. Finally, we apply the corrected statistics to independence 
tests on epidemiologic and ecological data. 

1 Introduction and notations 

Physical, sociological or biological surveys often lead to data presented as 
contingency tables. These surveys aim at studying relationships such as 
total, mutual, partial or conditional independence between several characters 
in a population. Table notations are very useful to formulate the test and 
make the corresponding hypotheses explicit, but can be cumbersome when 
the number of characters exceeds three. For theoretical results, we therefore 
use vector notations instead, and we reformulate the independence test as a 
multinomial goodness of fit test. 



1.1 Goodness of fit tests 

Let p = (pi, . . . ,Pr) be a probability vector of dimension R where R ^ 2 
is the total number of cross-classifying categories. Let n be the sample 
size and x = (rti, . . . ,nn) the vector of observed frequencies, realization of 
X = (N\, . . . ,Nr) with distribution M(n;p), multinomial distribution with 
parameters n and p. We denote by p° the probability vector under the null 



1 



hypothesis and consider the following test 



Ho- p = p° against Hi : p / p° . (1.1) 

Popular goodness of fit statistics include Pearson's chi-square statistic 
Q and Kullback's minimum discrimination information statistic G, defined 
in [10]. For two probability distributions p and p' , these are written : 

Q P (p')=nJ2 iP ' r ~ Pr)2 , (1.2) 
P- 

and 

G p (p') = 2nJ2p>-- (1-3) 

They belong to the power divergence statistics family {i?C A , A £ M} defined 
by Read and Cressie in [5], respectively for A = 1 and A — > 0, where : 



9n _ 
r=l 



/ \ A 
Pr 



For a review on goodness of fit testing methods and statistics, see [3]. 

The vector p° is not always completely specified. We assume that p° is 
a known function of an unknown parameter 9 of C M s with s < R — 1, 
and denote p° = p (0) with = (0i, . . . ,9 S ). We suppose that the functions 
9 i— > we consider here are bijective and estimate p°(9) by p*° = p°(9*) 

where 9* is the maximum likelihood estimator of 9. As the probability p 
is generally unknown, we also estimate the p r for r in {1, . . . , R} by their 
maximum likelihood estimators p* = n r /n. Throughout this paper, under- 
lying indexes n are omitted for simplicity of notation, and we consider the 
statistics Q p *o(p*) and G p *o(p*). 

Read and Cressie show that under Birch's regularity conditions, see [2], 
the statistics RC X are asymptotically equivalent and share a common chi- 
square limit distribution: 

Theoreme 1. 

VAel, lim £ no (RCL(p*)) = Xr-s-i- (1-5) 

n— >+oo ^ 

This implies the following result : 
Corollaire 1. 

lim £ho(Q p *°(p*)) = liin £h {G p *°(p*)) = Xr-s-i- ( L6 ) 
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This asymptotic result is a consequence of the Central Limit Theorem. 
Classical empirical limitations include that the sample size n must be over 30 
and that all expected frequencies must be over 5. Usually G p *o(p*) is pre- 
ferred to Q p *o(p*) because it is less sensitive to small cell frequencies. Read 

and Cressie recommend the use of RC p i (p*) instead for n 10 and a min- 
imum expected frequency over 1. Sometimes, even the lower bound 0.5 is 
accepted for the expected frequencies. A review on these conditions can be 
found in 

1.2 Sparse tables 

When there are too few subjects in the study or when the classifying cate- 
gories are too numerous, the table comprises one or several empty cells called 
random zeros. The table is then called sparse and it is likely that at least one 
cell has an expected frequency below 0.5. Random zeros would not appear 
if the sample was of sufficient size. Structural zeros, corresponding to cells 
with an expected probability of zero, are not considered here and should be 
suppressed. We therefore assume that the following condition is satisfied : 

pf^O, r €{1,...,R}. (1.7) 

Sparse tables can nonetheless be tested for independence, for example by 
regrouping cells so that the condition on the expected frequencies is satis- 
fied. However this procedure is not always data relevant. Fisher's exact test 
given in [6] applies without restrictions, except that it becomes numerically 
unmanageable when the table dimension grows. This leads to the use of 
Monte Carlo simulation methods as explained in [T]. 

The approach we propose here consists in correcting the historical statis- 
tics Q p *o(p*) and G p ,o{p*) according to the number of zero cells, by gener- 
alizing and improving a method designed by Ku in [9]. 

2 Corrections for Pearson's and Kullback's statis- 
tics 

Let C be the random variable giving the number of zeros in the vector 
or the contingency table, and c a realization of C with R — c ^ 1. For 
simplicity, we assume that ri\ = n% = ■ ■ ■ = n c = and that nj 
1 for all j in {c + 1, ...,12}. The maximum likelihood estimator p* = 
(0, . . . , 0, n c+ i/n, . . . , nn/ri) underestimates the pi for i in {1, . . . , c} and 
overestimates the pj for j in {c + 1, . . . , R}. Its use when c ^ thus has 
consequences on the statistics Q p *o(p*) and G p *o(p*). 
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2.1 Ku's correction for one zero 

Ku argues in [9] that G p *o{p*) tends to innate with respect to Q p *o(p*) when 
C grows. He then proposes to subtract 1 from G p *o(p*) for each zero, that 
is c in total. He proves the asymptotic equivalence between G p *o(p*) and 
its corrected version only for c = 1. I will explain why his reasoning is 
inconsistent. First, he considers a new statistic which is not a member of 
the power divergence family. Moreover, he uses the straightforward Lemma 
[I] to deduce the approximation: 

\ np* u J np* v 

for a = n r /n, b = p*° and n r = 0, despite the fact that a = and 1/a is not 
bounded. 

Lemme 1. For each a,b > such that a < 2b, b < 2a and the quantities 
1/a and 1/6 are bounded, we have : 

. (a\ a 2 -b 2 , 

The expression on the left in (|2.ip is null whereas the one on the right 
is negative. The sum over r of the left-hand side gives G p *o{p*) and we 
recognize Q p *o(p*) in the sum of the right-hand side. However, unlike Ku 
seems to think, zero and non-zero cells tend to compensate for each other 
and the approximation (|2.ip can not be extended to the corresponding sum. 
There is therefore no behavior of G and Q that we can deduce from this. 
Finally, he illustrates this asymptotic result on a small sample of size n = 10. 

We propose new corrections for both statistics, based on Ku's correction 
and a likelihood inequality that we present in the next subsection. 



2.2 Likelihood inequality 

Let us consider the following inequality coming from a likelihood reasoning : 
the sample vector we observe can be thought of more likely to happen than 
any other possible vector, since it is the one we actually observed. With 
exactly c zeros and R — c non-zero cells observed, so for all m ^ c and n'- 

such that n'j < nj, Vj £ {c+1, ... , R} and n = Y^=l n 'i + Sf=c+i n 'j '■ 

P(Aq = 0,...,JV c = 0,iV c +i = n c+1 , ...,N R = n R )> 
P(JVi = n' 1 ,...,N m = n' m , N m+1 = . . . = N c = 0, iV c+ i = n ' c+1 , ...,N R = n' R ). 

(2.2) 

We give in Proposition [1] a sufficient condition on the p r for (|2.2|) to be 
satisfied, and prove this statement in Appendix [A] 
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Proposition 1. The inequality (|2,2p is satisfied under the following assump- 
tion : 

Pi^^ 1 , Vie {l,...,c}, Vj€{c+l,...,i2}. (2.3) 
n 

2.3 Corrections 

We propose an estimator p for p, different from the maximum likelihood 
estimator, with the following form : 

(Pi =a, Vi G {l,...,c}, 

\pj =^;-d, Vje{c+1,...,R}, 

where a = a n , b = b n and d = d n are random variables depending on n 
designed to compensate for the under- and overestimations due to p*. We 
thus take them positive with < b < 1, such that b inflates the modified 
maximum likelihood estimator nj/n b and such that d controls the related 
rise. 

This new probability vector allows us to define corrected statistics Q p *o(p) 
and G p *o(p), provided that the parameters a, b and d satisfy several con- 
ditions. Forcing the summation of the p r to 1 implies that (R — c)d = 
ac + n l ~ b — 1, and thus allows us to define Q p *o(p ab ) and G p ,o(p ab ) with p ab 
such that : 



Definition 1. 



ac + n 



1-6 



Vi G {1, . . . ,c}, 

1 (2-5) 
-, Vj e{c + l,...,R}. 



n b R-c 
Note that for c = 0, we fix a = and 6=1. 

Before considering the other conditions on p ab , let us first give some 
notations. Let n be min Jg .r c _ t _i i . ..,_??}{%} an d n be maxj g { c+lj ^{n.,}. Let 
n stand for n — n(R — c) and n for n{R — c) — n. We dismiss the uniformly 
distributed case where : 

IX 

n = n = nj = — , V? G {c + 1, . . . , R}, (2.6) 

it — c 

thus guaranteeing that n and n are positive. Let 

/ ln(fl/(fl-l)) ln(n) ln(n-n) ^ 
6 min = max 0, --— ,T~T^^ i i \ — ' ( 2 - 7 ) 

( (n-n b )(R-c)+n b \ , . 

a min {b) = max 0, ^ A ; , 2.8 

V en" / 



5 



and 



n b — n n b — n 



« max (6) = min 1, —5= _ - - - , (2.9) 



Proposition [T] is applied to p ab . Together with inequalities < < 1, for 
all r in {1, ... , R} it is then equivalent to these conditions on a and b: 

b m in <b< 6 max = 1, (2.10) 

and 

^min 

(b) < a < a max (6). (2.11) 

We want to make a practical choice among possible values of p ab ensuring 
us that it is as far from p* as possible. We therefore fix b in (|2.10p quite far 
from 1 as a convex combination of 6 m i n and 6 max with an empirical parameter 
h equal to 0.1. For this value of b we choose a near the upper limit of the 
interval in (|2.1ip . that is : 

b = /ifemax + (1 - h)b min and a = a max (b) - e, (2.12) 

where e is a small constant designed to eliminate boundary effects. 

The final expressions we get for the corrected statistics Q ab = Q p *o(p ab ) 
and G ab = G p *o{p ab ) of Q = Q p *o{p*) and G = G p *o(p*) are: 

Q ab = n 2 ^Q- f(a,b), (2.13) 

with 



/(a, 6)= n(l-n^) + 2 " 1 -V^->-l) g 



-a 



R — c ^— ' np*; 

j=c+l J 
c 1 / , 1-6 1 \ 2 R 

2 1 I ac + n — 1 N 



and 
where 



G ab = n 1 ~ b G - g(a,b), (2.15) 



g(aM _ 2n | _____ ^____i__L_ 1 



a Vlnf— ^ -n 1 - 6 V Uj In f - c) - + — - - ^ ] 



(2.16) 



2.4 Convergence 

In this paragraph, we show the convergence in distribution of Q ab and G ab 
to a chi-square distribution. Let us first study the parameter b. 

Propriety 1. Bound b m i n is strictly less than 1. Moreover, if n = o(n) and 
n ~ n when n tends to +oo ; then b m i n and b tend to 1. 

Proof. Inequalities n{R — c) < nR, n — n(R — c) < n and n — n<n show that 
all three components of the maximum defining 6 m j n are strictly less than 1. 
Their order 1 developments as n tends to +oo give their convergence to 1, 
with the quantities R — 1, R — c and R — c — 1 bounded. □ 

We also study the asymptotic behavior of C, for once denoted C n , and 
state the following lemma which proof appears in Appendix lAl 

Lemme 2. The number of zeros C n converges almost surely to as n tends 
to +oo. 

We recalled in section [1] the convergence of Pearson's and Kullback's 
statistics to a chi-square distribution. In Theorem [2j we state a similar 
property for the corrected statistics Q ab and G ab . 

Theoreme 2. Under Birch's regularity conditions in f^, the estimator p ab 
defined by , (f!H0|) and f|27TT]) is such that : 

lim C Ho (Q ab )= lim C no (G ab ) = x 2 R -s-i- (2-17) 

n— >+oo n— >+oo 

Proof. We deduce from Lemma [2] the existence of a set Q' of probability 1 
on which we can find a rank no such that C n = for n ^ hq. The variable a 
is then set equal to and the variable b equal to 1. Hence, the estimates p* 
and p ab match, so that Q ab = Q and G ab = G. Theorem [1] then completes 
the proof. □ 

The final two sections are dedicated to proving the relevancy of our cor- 
rections through simulations and real data analyses. 



3 Simulations 

In this section, simulations confirm the necessity to correct not only G but 
also Q. We compute the statistics Q, G, RC 2 ^ = RC 2/ J(p*), Q ab and G ab 
on 1000 vectors of length R = 100 of total frequency n = 400 for each of the 
four multinomial distributions defined in Table [T]by f\ to f^. Let E r denote 
the expected frequencies for r in {1, . . . , 100}. 

For each set of vectors sharing the same c, we give the quantile of order 
1 — a = 95% for the five statistics considered. Results are displayed in 
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Table 1: 


Multinomial probabilities fx to f±. 




|{r; E r < 0.5}| 


Probability 


h 


20 


(0.0002, . . . , 0.0002, 0.01245, . . . , 0.01245) 

N v ' V v ' 


h 


50 


20 times 80 times 

(0.0002, . . . , 0.0002, 0.0198, . . . , 0.0198) 


h 


70 


50 times 50 times 

(0.0002, . . . , 0.0002, 0.03286667, . . . , 0.03286667) 

V v ' " v ' 


h 


90 


70 times 30 times 

(0.0002, . . . , 0.0002, 0.0982, . . . , 0.0982) 

90 times 10 times 



figure [TJ as well as a line indicating the chi-square quantile of order 1 — a = 
95% with R— 1 = 99 degrees of freedom Xo 95 99 = 123.22. Only the center of 
each graph should be considered because the quantiles for extreme numbers 
of zeros are computed on too few observations, sometimes only on one or 
two among the 1000 simulations in total. 

Quantile values of Q tend to explode as c grows whereas G stays quite 
stable around ^09599- This behavior is the opposite of the one predicted 
by Ku. For fx, statistics Q ab and G ab lead to the rejection of the null 
hypothesis. For fa to f± however, their quantiles lie below the critical line 
and Hq is accepted. We thus have compensated for the rise of Q, and both 
our corrected statistics are stable. 

This analysis is confirmed by the computation of empirical risks of type I 
for a = 0.01, 0.05 and 0.1 as showed in Table [2j and by the power study 
below. Probabilities fx to f± are perturbed into f[ to f^ such that for all j 
in {1,2,3,4}: 

Vi€{l,... I 10}, f'j(i) = ^(0 + 1/300, (3.1) 
Vi€ {11,..., 90}, /j(z) = /,•(»), (3.2) 
ViG {91,..., 100}, /<(*) = 1/300. (3.3) 

Vectors are simulated with probabilities fx to f±, and goodness of fit for 
f{ t° fi is tested. The Tables [2] and [3] show that the empirical type I risks are 
lower for our corrections compared to the classical statistics when c is quite 
important, whereas the empirical power is much higher for our corrections 
when c is small. 

4 Applications 

We apply the total independence test using the corrected statistics Q ab and 
G ab to two datasets involving two-dimensional tables. For such tables, the 
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14 16 18 20 22 

C 








Figure 1: Quantiles of order 0.95 for Q, Q ab , G, G ah and RC 2 / 3 as functions 
of c, under null multinomial probabilities f\ to f^, for 1000 samples of size 
n = 400 and R = 100 categories. The line represents the threshold Xq 95 99- 
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Table 2: Empirical type I risks for Q, Q ab , G, G ab and RC 2 / 3 , for 1000 
samples of size n = 400 and R = 100 categories, at levels a = 0.01, 0.05, 
0.1, for vector probabilities f\ to f^. 



o 


n 


mode(c) 


Q 


Q ab 


G 


Qab 


RC 2 / 3 


0.01 


h 


19 


0.031 


0.233 


0.003 


0.401 


0.003 




f-2 


47 


0.030 


0.005 













h 


65 


0.027 
















h 


84 


0.031 














0.05 


h 


19 


0.044 


0.352 


0.010 


0.573 


0.017 




h 


47 


0.024 


0.005 













h 


65 


0.069 











0.006 




fi 


84 


0.052 














0.1 


h 


19 


0.130 


0.479 


0.033 


0.674 


0.039 




h 


47 


0.072 


0.005 





0.005 







h 


65 


0.137 
















fi 


84 


0.098 











0.006 



Table 3: Empirical powers for Q, Q ab , G, G ab and RC 2 / 3 , for 1000 samples 
of size n = 400 and R = 100 categories, at levels a = 0.05 for simulated 
vector probabilities fx to and null vector probabilities f[ to f^. 



U 


Hi 


mode(c) 


Q 


Q ab 


G 


Qab 


RC 2 ' 3 


fi 


fi 


19 


0.233 


0.853 


0.322 


0.983 


0.157 


h 


fi 


47 


0.087 


0.026 


0.009 


0.061 


0.009 


h 


fz 


64 


0.229 


0.005 








0.027 


h 


fi 


84 


0.089 
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Table 4: Diplotype table for the association between TNFAIP3 and 

Systemic Sclerosis.. 



Status HI/HI H1/H2 H1/H3 H1/H4 H1/H5 H1/H6 

Sound 98 7 116 2 71 3 

Affected 91 9 104 3 70 12 

H2/H3 H2/H5 H2/H6 H3/H3 H3/H4 H3/H5 

Sound 4 2 34 I 42 

Affected 5 4 1 30 2 40 

H3/H6 H4/H5 H5/H5 H5/H6 

Sound 2 1 13 1 

Affected 7 1 13 5 



hypotheses are written : 

%o : Pij=Pi+P+j against Hi : 3 (ii, ji) € I x J, ^ p h+ p +jl , 

(4.1) 

where £>j + and p + j are the marginal distributions for the two characters 
featured in the table. To ensure the condition (|1.7p we remove the empty 
lines i of {1, ...,/} such that = 0, and the empty columns j of {1, ... , J} 
such that n+j = 0. 

4.1 Multi-marker approach for Systemic Sclerosis 

Table |4] is the diplotypic table obtained from an association study in Humans 
looking for an association between three genetic markers on the gene TN- 
FAIP3 and Systemic Sclerosis presented in |S]. Empty columns have been 
removed. A haplotype is the allelic distribution of markers on a chromosome, 
and a diplotype is the combination of both parental haplotypes. Diplotypic 
tables, though more interesting than haplotypic tables because they take 
into account more information, are usually trickier to handle because they 
are sparse. Our corrected statistics can therefore be helpful in such situa- 
tions. 

Though nine haplotypes theoretically exist, only eight are observed, de- 
noted HI to H8, leading to 8 2 = 64 diplotypes Hi/Hj. Two samples are 
compared, affected versus sound subjects, on which we test the indepen- 
dence between the diplotype configuration and the health status of n = 794 
individuals. 

The table is of dimension 2 x 16, that is R = 32 categories with c = 1 zero 
and s = 16 parameters. There are exactly 16 expected frequencies below 5. 
The chi-square quantile Xo 95 15 = 24.99 is compared to the statistics : 

Q = 14.62, Q ab = 20.76, G = 15.82, G ab = 28.43, RC 2/3 = 14.85. 
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Table 5: Contingency table for the joint study of trophic level and 
vegetable composition in rivers. 



Trophic level 


(0,0,0) 


(1,0,0) 


(r,P-, 
(0,1,0) 


<0 

(0,0,1) 


(1,1,0) 


(0,1,1) 


Oligotrophic 








3 





3 


2 


Mesotrophic 


2 


1 





2 


1 





Eutrophic 


2 





3 


1 


1 






Only G leads to reject the null hypothesis of independence. This seems to 
be the right decision since it is confirmed by the single-markers approaches 
and haplotype tests in |5], all showing a significative association between the 
markers and the disease. 

4.2 Trophic level and vegetables in the rivers of the Petite 
Camargue Alsacienne 

The search for a link between the trophic level and the vegetable composition 
of some rivers of the Petite Camargue Alsacienne in North-East France leads 
to Table [5] A river can be either oligotrophic, mesotrophic or eutrophic, if 
its nutritive content is respectively poor, intermediate or high. Uncommon 
vegetables are considered rare, exotic or polluto-tolerant. To each river, a 
triplet of binary characteristics (r, p, e) is assigned, indicating the presence 
(1) or the absence (0) of rare (r), exotic (e) and polluto-tolerant (p) species. 
The original ecological study can be found in |11| . We consider n = 21 
different rivers. Two empty columns were removed from the original table, 
leading to Table [5] of dimension 3x6, that is R = 18 categories and c = 7 
zeros, with s = 17 parameters. Here are 3 expected frequencies below 0.5. 
Test statistics are compared to the chi-square quantile Xo 95 10 = 18.31: 

Q = 14.38, Q ab = 20.68, G = 18.67, G ab = 26.05, RC 2/3 = 14.84. 

Both corrected statistics Q ab and G ab as well as G lead to reject the null 
hypothesis, indicating an association between trophic level and vegetable 
composition. A thorough study of the table shows that rare species tend to 
settle preferentially in oligotrophic rivers. They are indeed better adapted 
to this kind of environment which tends to disappear from the rivers in 
the study. Moreover, polluto-tolerant species constitute the majority of the 
vegetables in eutrophic rivers. A eutrophic environment is competitive and 
these resistant species tend to get the best of it. 

4.3 Discussion 

We suggest to compute both Q ab and G ab , and to reject the null hypothesis 
if at least one of them is larger than the chi-square threshold. 
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Our results tend to prove that this approach is relevant and our correc- 
tions efficient in sparse tables. They are all the more interesting for the fact 
that sparse tables are usually left aside because the hypotheses needed to 
apply classical chi-square tests are not satisfied. 



A Appendix section 

Proof of Proposition^ Assume that the first m of the c frequencies n,;, 1 ^ 
i ^ c are modified. Let Tipj £ {c+1, ... , R} compensate for these changes. 
The likelihood inequality (|2.2p is then equivalent to : 

(n c+ i-n' c+1 ) (n R -n' ) n < n > 

Pc + i x x Pr > p i ... Pm . (A.l) 

(n c+ i - n' c+1 + 1)! (n R - n' R + 1)! ^ n[\ n' m \ 

Let us show that f j 2 . 3 1) implies (lA.ip . Applying (|2.3p to each element of the 

R m 

left hand side of fjA. 1[) with multiplicities such that : (rij —n'j) = n i: 

j=c+l i=l 

we get : 

(n c+1 -n' c+1 ) (nR-n'n) 



n (n c+1 -n' 0+1 +l) n {n R -n' R +\) 

As (rij -n'j + iy. ^ n {nj ~ n '^ for all j in {c + 1, . . . , R} and n<! ^ 1 for all i 
in {1, . . . , m}, we deduce fjA. 1[) and equivalently (|2.2p . □ 

Proof of Lemma\^ Let us first show that: 

Ve > 0, lim F(C n > e) = 0. (A.3) 

n— y +oo 

For each n we compute P(C n = c) for c in {0, . . . , R — 1}. Let p° be the 
probability under the null hypothesis p° = (p®, . . . ,p®, . . . ,P%)- A 

subject belongs to one of the first c cells with probability q° = p\ + • • • + p® 
and to one of the R — c last cells with probability 1 — q° = p°+i + ' ' ' + Pr> 
with q° in ]0, 1[. 

Let us consider now the binomial distribution B(n; q°) and write the 
probability P(C n = c) of obtaining a table containing exactly c zeros placed 
anywhere : 

F(C n = c) = (f\v(Ni = ... = N c = 0, N c+l ^0,..., N R ^,A) 
= Q(1-«°) B . (A.5) 

Then : 

P(C n = 0) = 1 - = c) = 1 - P(C n > e), Ve e]0, 1[. (A.6) 



c=l 
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As (^) is bounded and (1 — q°) n tends to 0, the probability P(C n = c) also 
converges to for all 1 ^ c ^ R — 1 when n tends to +oo, and so does the 
corresponding sum over c. 

We now use Borel-Cantelli's Lemma to conclude that C n converges to 
almost surely. Indeed, for e > 0: 

R— 1 

£ ncn > £ kx; nc n > i) = £ ( R c ) ^ < +°°- ( A - 7 ) 

nSsl n>l c=l / * 

□ 
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